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We construct energy-optimized resonating valence bond wavefunctions as a means to sketch out 
the zero-temperature phase diagram of the square-lattice quantum Heisenberg model with competing 
nearest- (Ji) and next-nearest-neighbour (J2) interactions. Our emphasis is not on achieving an 
accurate representation of the magnetically disordered intermediate phase (centred on a relative 
coupling g = J2/J1 ~ 1/2 and whose exact nature is still controversial) but on exploring whether 
and how the Marshall sign structure breaks down in the vicinity of the phase boundaries. Numerical 
evaluation of two- and four-spin correlation functions is carried out stochastically using a worm 
algorithm that has been modified to operate in either of two modes: one in which the sublattice 
labelling is fixed beforehand and another in which the worm manipulates the current labelling so as to 
sample various sign conventions. Our results suggest that the disordered phase evolves continuously 
out of the (71-, 7r) Neel phase and largely inherits its Marshall sign structure; on the other hand, the 
transition from the magnetically ordered (tv, 0) phase is strongly first order and involves an abrupt 
change in the sign structure and spatial symmetry as the results of a level crossing. 



I. INTRODUCTION 

Simple spin models have contributed significantly to 
our understanding of quantum magnetism. They consist 
of mutually interacting spin-S" objects arranged in a lat- 
tice and are meant to describe the behavior of localized 
electrons in a crystalline environment. Such models are 
generally viewed as effective, low-energy descriptions, de- 
scended from their electronic parent models by a process 
of integrating out the gapped charge degrees of freedom.^ 

A tremendous variety of spin interactions can arise. In 
particular, a "t/f/"-style powerseries from the strong cor- 
relation limit generates (or at least motivates) an increas- 
ingly complicated zoo of multi-spin interaction termsPHU 
Nonetheless, we know that even the leading order term 
in the expansion, corresponding to Heisenberg models 
with just two-spin interactions, can display highly non- 
trivial physics if the exchange interactions are sufficiently 
frustrating.*^ l n that case, the ground state may be a 
magnetically disorde red, spin-rotation-invariant state — 
either liquicPor soli d 10 * 11 -* — having no classical analogue. 

Otherwise, conventional magnetic order (at some or- 
dering vector Q) is a generic feature of the ground state 
for Heise nberg models in spatial dimension greater than 
one) 12 * 13 * The absence of frustration is connected to three 
inter-related properties: (i) the existence of a bipartite 
labelling such that all antifcrromagnetic interactions con- 
nect sites in opposite sublattices, (ii) strict adherence to 
a Marshall sign rulep* and (iii) the possibility of trans- 
forming mechanistically to a basis in which all amplitudes 
of the wavefunction are real and nonnegative. The last 
of these is why nonfrustrated models can be easily simu- 
lated using quantum Monte Carlo approaches ESHI3 

For the S — 1/2 case, all three properties are con- 
ceptually unified in the language of valence bondsW^^ 
The collinear, Q-ordered ground state of a nonfrustrated 



Heisenber g mod el can be described in a bipartite valence 
bond basi d 22 * 23 -* in which the AB sublattice labelling coin- 
cides with the alternating pattern laid out by Q and only 
spins in opposite sublattices are bound into singlet pairs. 
In terms of such a basis Vab = {\v)}> the ground state 
has an expansion |?/>) = ip(v)\v) in which each ampli- 
tude ij)(v) is real and nonnegative. The exact amplitudes 
can be obtained numerically by projectionl^HIII 

It is also possible to find extremely good approximate 
values of the form ip(v) ~ FJ^ /i(ry), where h{v) > 
is a function of the vector connecting bond end points. 
This resonating valence bond (RVB) ansatz, due to 
Liang, Doucot, and Anderson,^ strictly enforces the ge- 
ometric tiling constraint on the singlet bonds but ignores 
additional bond-bond correlations.^ For a magnetically 
ordered state, one can show that factorizability int o in- 
dividual bond amplitudes is the correct assumption! 30 * 31 1 
Moreover, for nonfrustrated systems, the amplitudes ex- 
hibit powerlaw decay, and hence the wavefunction con- 
tains bonds on all length scales. 

As a specific and illustrative example, we consider the 
square-lattice J1-J2 model for spin half. It has two 
nonfrustrated limits. The model with anitferromagnetic 
nearest-neighbour interactions only (Ji = \,Ji = 0) 
exhibits a Nccl ordered ground state whose staggered 
moment is roughly 60% of its fully polarized, classi- 
cal value. The state is almost perfectly captured by 
an RVB wavefunction whose bond amplitudes are com- 
puted as h(r) = £ q e lqr [l - (1 - 7q ) 1/2 ]/7 q - Here, 
7 q = (cosq x -\-cosq y )/2, and the wavevector sum is taken 
over a Brillouin zone reduced with respect to Q = (tt, tt). 
The opposite limit, with nea;i-nearest-neighbour interac- 
tions dominating ( J\ = + , J 2 = 1), is equivalent to two 
interpenetrating nearest-neighbour Heisenberg antiferro- 
magnets rotated 45° . The spin directions in the two oth- 
erwise disjoint subsystems lock to each other^l provided 
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that Ji is not strictly zero. In this case, the ground state 
is equally well described by the RVB wavefunction, but 
with the substitution of 7 q = cos q x cos q y and a Brillouin 
zone defined modulo Q = (n, 0) or Q = (0, it). 

What we present in this paper is an attempt to in- 
terpolate between these two limits — through the entire 
range of relative couplings that are highly frustrated — 
using the RVB state as a variational wavefunction. Our 
approach is inspired by Ref. [331 but there are several im- 
portant differences. The first is simply the scale of the 
calculation: we have simulated a large number of lattice 
sizes up to L = 32 on a dense grid of relative coupling 
values [g — J2/J1 ranging from to 1 in steps of 0.01). 
Second, we do not require that h(r) respect the full C4 
symmetry of the square lattice. Rather, we impose only 
the x- and y-axis reflection symmetry, giving the ampli- 
tudes an opportunity either to acquire (over the course 
of the energy optimization) the full symmetry or to set- 
tle into a state that looks different under 90° rotation. 
Third, we explore the space of AB sublattice labellings 
by which the bipartite valence bond basis is constructed. 

As in Ref. [33J we make use of an unbiased, stochas- 
tic optimization scheme. Changes to the h(r) values are 
made according to the sign of the local energy gradient. 
Step sizes are randomized, and their magnitude decreases 
on a powerlaw schedule. We do not attempt to guide 
the optimization, other than to ensure that none of the 
bond amplitudes goes negative; nor do we impose any 
constraints on the variational parameters based on a ny 
prior knowledge (gleaned, e.g., from mean field theory 3 -^ 
or from a master equation analysis^). 

We discover the following. At this level of approxima- 
tion, the J\-Ji model does indeed support a magnetically 
disordered intermediate phase. But its width is much 
smaller than expected: the phase boundaries are found 
to be at g c i = 0.54(1) and g c2 = 0.5891(3). The transi- 
tions are unambiguously second- and first-order, respec- 
tively, with the ground state always achieving the full C4 
symmetry for g < g c i. As the system is tuned up from 
g = 0, increasing frustration eventually extinguishes the 
(tt, 7r) ordered moment at g c \ in a continuous fashion. 

The disappearance of magnetic order is preceded by 
a failure of the Marshall sign rule at <?mi = 0.398(4). 
Still, even though the rule is not strictly obeyed beyond 
<7mi, the Marshall structure inherited from the g = 
model remains largely in tact throughout the intermedi- 
ate phase. This is true in the sense that continuing to de- 
fine the bipartite bond basis from a checkerboard sublat- 
tice decomposition produces only a microscopic number 
of negative h(r) values — only /i(±l,±2) and /i(±2,±l) 
initially. Moreover, when we allow the AB pattern to 
arise on its own within the simulation (described in de- 
tail in Secs.|IID and III), the checkerboard pattern is the 



one selected whenever g < g C 2. 

On the other hand, the RVB state at large g explicitly 
breaks the 90° rotation symmetry and has a Marshall 
sign structure based on a stripe sublattice decomposi- 
tion. As the coupling is tuned down from the g = 00 



limit, the (7r, 0) ordered moment is not strongly affected, 
and it persists with only weak variation (never dropping 
below 47% of its fully polarized value) down to g c2 , where 
the spatially symmetric, checkerboard-based RVB wave- 
function takes over as the lowest energy state. This state 
in the region g c \ < g < g c i is, as far as we can tell, 
featureless. It exhibits no long-range spin or dimer or- 
der, and it breaks no symmetries. It is not, however, a 
"short range RVB state" in the usual sense, since it is 
not made up of predominantly short bonds. Its ampli- 
tude function h(r) is highly anisotropic (as anticipated 
elsewhere^) and remains long ranged along the princi- 
ple spatial axes. Spin correlations appear to be criti- 
cal, and dimer correlations decay either exponentially or 
with a high power law. This is in stark contrast to the 
usual short-bond-only RVB state, which has spin corre- 
lations that decay expo nentia lly^ and dimer correlations 
that decay algebraically! 35 * 36 ^ The presence of long bonds 
implies an absence of the topological ordeiP^S that is 
characteristic of the purely short ranged RVB state in 
two dimensions. 



II. MODEL AND METHOD 

A. Frustrated Hamiltonian 

The spin-half, square-lattice Heisenberg model with 
frustrating interactions has a Hamiltonian 



H = Ji Si 



{{hi)) 



(1) 



where J\ > and J2 > are the antiferromagnetic ex- 
change couplings. The summations range over pairs of 
adjacent sites and over farther pairs that sit 

diagonally across a plaquette. The ratio g — J2/J1 is 
the key tuning parameter at zero temperature. In the 
classical version of this model (5 — > 00), two magnetic 

phases meet at exactly g — 0.5, separated by a first-order 
transitionPHSD 

In the 5=1/2 problem, the two magnetically ordered 
ground states obtain for values g < 0.4 and g > 0.6p3El 
and a magnetically disordered phase intervenes. (There 
is, however, a good deal of disagreement over the exact 
positions of the critical points; cf. Refs. 1431 and 1551 which 
put the lower critical point as low as 0.35 and as high as 
0.45.) The physics of the phase in the intermediate region 
is not known with complete certainty, but it is believed 
to be short ranged and not to exhibit any kind of con- 
ventional magnetic order. One possibility is a crystalline 
arrangement of valence bonds, a state with broken trans- 
lational symmetry in which singlet formation favours an 
enlargement of the unit cell beyond that of the underlying 
square latticeP^HZHI A featureless spin liquid that doesn't 
break any s ymme tries was also proposed as a possible 
ground state.™™ 

The case for a spin liquid has been advanced by re- 
cent tensor produclP^ and density matrix renormaliza- 



3 



O O O O O'" 

OiiiiiiQiiiiiiQiiiiiiQiiiiiiQiu 



O 0"\"0 O"' 

QIIIIIIQIIIIIIQIIIIIIQIIIIIIQIII 



O Cf^O 6 Q" 

QiiiiiiQiiiiiiQiiiiiiQiiiiiiQiii 



O O O O O'" 



o o 



o o 

6 6 



o o 

§ <5 

o 6 

6 6 



o p„ o a o p, o 
\ # \ v % 
q, o n o p o o 



o 3^-0 X °M ° 

% \ />- \ 

oVoVoVo 



FIG. 1. (color online) Dashed lines indicate the Ji (red) and 
J2 (blue) exchange couplings. The basis contains only prod- 
uct states of singlets connecting sites in opposite sublattices. 
(Left) In the limit g — J2/J1 = 0, a checkerboard pattern 
of A and B labels that coincides with (n, n) magnetic order. 
(Right) In the limit g — 00, a stripe pattern that coincides 
with (ir, 0) order. In each case, three permissible singlet pair- 
ings are indicated. 



tion group calculations.^ With regard to the latter re- 
sult, Sandvik has suggested that the use of a cylindrical 
geometry complicates the detection of crystalline order 
His numerical experiments seem to indicate that the mix- 
ture of open and closed boundary conditions significantly 
raises the crossover length scale £ beyond which bond or- 
der takes hold (i.e., where the finite size scaling behaviour 
of the dimer-dimer correlations is truly in the asymptotic 
regime). Such questions are difficult to resolve. Unlike 
in three-dimensional systems, where crystalline bond or- 
der order, if it is present, is almost always strongP^ 
in two dimensions it is quite delicate and can easily be 
disguised by a U(l) effective symmetry for system sizes 
L ~ £• (See Sects. Ill and IV of Ref. (SO! and references 
therein.) Here, we attempt to make the best of this un- 
satisfactory state of affairs. We simply take the point of 
view that, for the lattice sizes (up to L = 32) we can 
simulate, the liquid and the weakly ordered bond crystal 
are indistinguishable. 



B. RVB trial wavefunction 

In quantum Heisenberg models, competing interac- 
tions that frustrate the order have the potential to sta- 
bilize exotic quantum phases, but they also render the 
problem computationally intractable on large lattices. 
Frustrating interactions of even infinitesimal strength 
cause a sign problem — that makes quantum Monte Carlo 
calculations unfeasible. Moreover, the size of the Hilbert 
space grows exponentially with system size and is thus 
beyond the capability of exact diagonalization (ED) cal- 
culations if we want to get near the thermodynamic limit. 
(The record for spin half has recently jumped from 42 
sites™] to 48 sitesj^ a terribly impressive technical 
feat that nonetheless limits us to two-dimensional length 
scales ~ y48 that are quite small.) An approximate 



method based on good trial wavefunctions is therefore 
one of the few remaining possibilities for large systems. 

We consider a lattice of 27V spins and a factorizable 
RVB wavefunction of the form 



e n < 

v [i,j]£v 



(2) 



where the sum is over all partitions of the lattice into N 
directed pairs v = ([hJi], [12,32], [inJn])- To every 
such dimer covering v, there is a corresponding singlet 
product state; e.g., 
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in the 5=1/2 case. The set V = {|i>}} of all possible sin- 
glet product states is both overcomplete and nonorthog- 
onal and constitutes the so-called valence bond basis. 

We can now break up the lattice into two sublattices — 
groups of sites labelled A and B, equal in number — and 
restrict ourselves to a reduced basis in which valence 
bonds only connect sites in opposite sublattices (i.e., 
v £ V AB ^ S N , rather than v € V ~ S 2N /Z^). We 
adopt the convention that each bond [i,j] is arranged 
with site i in sublattice A and site j in sublattice B. 
This has the advantage of rendering the overlap strictly 
positive: (v\v') = 2 A ' ! ( C )~ JV , where Ni(C) is the number 
of loops in the double dimer covering C = {v,v'). (In 
this "bosonic" convention, the singlets are AB directed 
bonds. In the complementary "fermionic" convention, 
the bonds are directionless and all signs are moved into 
the overlapsPSiH) 

To start, we consider two families of trial state, each 
built using a bipartite bond basis consistent with one of 
two static choices of sublattice labelling, viz., the checker- 
board and stripe patterns shown in Fig. [l] Later in the 
paper, we go on to describe a procedure in which the 
trial state is built using an unrestricted bond basis and 
in which the sublattice labelling (and hence the Marshall 
sign convention) is determined dynamically. 

The RVB wavefunction is quite expressive. Its degrees 
of freedom are the full set of h(r) values with the bond 
vector r spanning all lengths and orientations that can 
be achieved on an L x L cluster with periodic bound- 
ary conditions and that are unique up to whatever sym- 
metries are enforced. (Still, the total number of pa- 
rameters grows only linearly with the number of spins, 
which is radically slower than the number of states in 
the total spin singlet sector.) Previous calculations of 
this kin d 33 * 34 ^ considered just the checkboard AB pat- 
tern and imposed on h(r) = h(x, y) the full symmetry 
of the lattice, such that h(x, y) = h(\x\, \y\) = h(\y\, \x\). 
In this calculation, we impose a less restrictive condi- 
tion, h(x,y) = h{[x\, \y\), that respects reflection sym- 
metry across the lines x — and y — but not across 
the lines y = ±x. For the checkerboard pattern, the 
number of free parameters is (L/2 — rf){L/2 — 1), where 
rj = (L/2 mod 2) distinguishes between L/2 even and 
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odd. For the stripe pattern, the count is only slightly 
higher: (L/2 + l)(L/2 - 1) = L 2 /4 - 1. 

To recapitulate, our work involves a basis choice. We 
do not construct the trial wave functions from the largest 
possible set of valence bond states in which the spins are 
joined in all possible ways. Instead, we obtain a more 
restricted basis by dividing the system into two groups 
of sites (A and B) and only keeping states in which bonds 
connect A sites and B sites (bipartite bonds) . No approx- 
imation is involved in this basis choice since the restricted 
basis is so massively overcomplete that even this subset 
still spans the relevant part of the Hilbert space. 

But in assigning A and B labels to the sites, we are 
making choices about the form of the trial wavefunction. 
By working with the checkerboard and stripe AB pat- 
terns, we are in essence adapting the trial wavefunction 
to g = and g = oo and taking advantage of the Mar- 
shall's sign rules that exist in the two limits. However, we 
are not biasing the wavefunction insofar as we are build- 
ing in any kind of magnetic order. The wavefunctions 
constructed from either AB pattern are fully capable of 
representing non-magnetic states. 



C. Sampling algorithm 

Every measurement (ip\0\ip) / (ip\ip) can be thought 
of as the ensemble average of an appropriate estimator 
O : C —> 0(C) in the gas of fluctuating loops described 
by 



1 

7n 



e^ i(c) n h > 



(4) 



Here, C = (v, v') is a loop configuration arising from the 
superposition of two dimer coverings, Ni(C) counts the 
number of loops, and q — 2 is the loop fugacity appro- 
priate for S = 1/2. When Marshall's theorem holds, the 
bond amplitudes satisfy hij > and thus every term in 
Eq. @ is nonnegative. This model is amenable to Monte 
Carlo simulation. We now outline a simple and efficient 
algorithm for performing the stochastic sampling. 

As a formal trick (in the spirit of Ref . , we enlarge 
the phase space from $o to $o x $i x ■ ■ ■ x where 
<I> n is the set of configurations in which 2n free endpoints 
have been introduced by breaking n valence bonds. (The 
system has been converted to one of both closed loops 
and open strings.) We take the partition function to be 



Z 



1 



n < 

[<j]ec 



(5) 



The configurations C are now assembled from all possi- 
ble partial coverings v = ([h,ji], [ia, 32], • ■ - > [injn]) of 
variable length < n < N, and a is introduced as a fu- 
gacity for the open strings (numbering N s (C) = N — n). 
The loop-only sector corresponds to the original parti- 
tion function, Zq = (!)$„. (In each string sector there 



is a Green's function defined by the string endpoints: 

Gij — (5i t a 1 Sj^fl 1 )c!> 1 , Gij-kl = (Si,a 1 Sj^j3 1 Sk,a 2 ^l,02)'^2t 

etc., and Z = (1)$ . Here a n and f3 n denote the po- 
sitions of the head and tail of the n th string. It is worth 
emphasizing that these 2n-point Green's functions do not 
coincide with expectation values of the physical spin op- 
erators. In general, we must take all measurements in the 
$0 configuration space using the loop esitmators derived 
in Ref. 122) 

We will consider a process that involves breaking a sin- 
gle valence bond ($0 — > $1) to produce an open string 
whose two endpoinds (the "head" and "tail") serve as 
walkers subject to Monte Carlo updates. The walkers 
move via a series of two-step motions that involve draw- 
ing a new bond and erasing an old one. When the walkers 
meet the loop is closed ($1 — > $o)- Figure [2] shows an ex- 
ample circuit. The fives successive steps shown in panels 
(b)-(f ) produce an overall change in the relative weight 
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Since we have chosen the bond amplitudes hij to 
be nonnegative, we can define a local amplitude Hi = 
J2j and a total overall amplitude H = J2i Hi = 
Yljj h^. These definitions will be useful in the deriva- 
tions that follow. 

To begin let us consider processes that take the system 
from the space of loops to the space of loops and one 
string. We move from a configuration C ~ \i, j] to a 
configuration C ~ (i)(J) by breaking a bond and 
thus leaving string endpoints (i) and (j). The transition 
probabilities for breaking and repairing the bond obey 
the detailed balance equation 



break , 



P(i)vc = W™f™P{3\i)*c> 



(7) 



Here P(i) is the probability of choosing a site i whose 
bond we want to break, and P(J\i) is the probability 
of choosing j given a walker (string endpoint) at site j. 
7r<7 and ttc 1 represent the likelihood of the system being 
found in configurations C and C . Their ratio is given by 



7Tc 



(8) 



If we choose which bond to break according to the dis- 
tribution of local bond weight P(i) = Hi/H and choose 
walker movements according to the distribution P(j\i) = 
/Hi, then 



TJ/brcak 
Tl/rcpair 



P{j\i)^c 

P{i)TT C 



a 
H 



(9) 



break 



We are free to choose a = H, in which case W, 
WfZff = 1. 

For motion of the walkers within $j, we need to know 
the transition rates between configurations C ~ k] 
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FIG. 2. (a) Two superimposed valence bond configurations form a collection of closed loops. Eight of the sites are numbered 
for use in Eq. Q. (b) Breaking one bond leaves an open string with a head and tail located at the former bond's endpoints. 
(c)-(e) The head and tail move by drawing a new bond and erasing the preexisting bond emerging from the destination site, 
(f) The open string is closed when the head and tail reconnect, (g) The repaired loop configuration, (h) Exchanging the 
background and foreground links in any loop is also a valid update. 



and C ~ [i,j](k). This represents a process in which 
a walker at i draws a new bond to some site j and then 
erases the prexisting bond connecting j to fc, thus leaving 
the walker at site k. The detailed balance equation is 



W?_XfP{j\i)Trc = W:^PU\k)n c 



The ratio 



*0 _ q SN, Kj_ 



(10) 



(11) 



depends on SN t = Ni(C')-Ni(C) = ±1 (or if the moves 
do not respect a fixed lattice bipartition; see discussion 
in Sect. II D I. As before, we attempt moves according to 



the distribution P(j\i) = hij/Hi. Then 



cwalk 



W{i -> fc) = P(J\k) ir c , = Hi 

W(k^l) P{j\i) TTC H k 



+ i 



(12) 



6/(1 



6) or 



which can be solved as WJ^t 
min(l, 5). 

Note that the transition rate does not depend on the 
ratio of bond amplitudes, as it would if we had, for ex- 
ample, selected a site uniformly with P(j\i) = 1/N. The 
ratio hij/hjk may fluctuate wildly over many orders of 
magnitude, so subsuming it into the sampling maximizes 
the efficiency of the algorithm. 

In the case of a translationally invariant system, the 
amplitude for pairing spins at i and j must be a function 
of the vector ry connecting the two sites; i.e., hij — 
ft(r y ). Hence, H = H/N = H t = £ r /i(r) for all i, 



which implies that P(i) — Hi/H — > 1/N is uniform and 
P(j\i) = hij/Hi — > h(rij)/H. The algorithm can be 
summarized as follows: 

1. Pick any valence bond (by choosing i uni- 
formly from the set of A sublattice sites and then 
selecting its partner site in v or v') and break it. 
The resulting string has endpoints at R = and 
R' = v 3 . 

2. To move the head, choose a new bond vector r from 
the distribution h(r)/H. So long as R+r ^ R', 
attempt to draw a new bond from R to R + r = 
(for some fc). The bond that already exists at that 
site [fc, I] is then erased and the walker is moved to 
r;. The move is accepted with probability 1 if it 
splits two loops and with probability 1/2 if it joins 
two loops. 

3. Otherwise, if R + r = R', close the open string by 
drawing a new valence bond from R to R'. 

The worm algorithm described here is ergodic and 
guaranteed to have a high acceptance rate. This is in 
contrast to the original bond-swap scheme proposed in 
Ref. [28l wherein two A-site or B-site bond endpoints 
sitting diagonally across a plaquette are swapped using 
Metropolis sampling. This antiquated algorithm runs 
into difficulty when the function h(r) is short ranged. In 
particular, short bonds that are adjacent but not shar- 
ing a common plaquette generate long bonds under rear- 
rangement, so whenever the amplitudes for long bonds 
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FIG. 3. (color online) (a) Two possible paths, marked 1 and 2, take the worm head to a site in the opposite or same sublattice 
of another closed loop. In either case, the loop is absorbed, (b) Path 1 leads to a rearrangement of the worm that preserves 
the AB labelling, (c) Path 2 requires that the AB labelling be reversed in the highlighted region, (d) Another worm, following 
two possible paths marked 3 and 4. (e) For path 3, the AB labelling is preserved, and the worm emits a new closed loop. (f,g) 
Path 4 requires that the AB labelling be reversed in the highlighted region. The number of loops remains unchanged. 



become small, the acceptance rate can become corre- 
spondingly small. Hence, ergodicity breaks down. Worse, 
there are typically many trapping configurations from 
which the simulation cannot emerge. The worm algo- 
rithm does not suffer from these problem, because it can 
traverse any local barriers by stepping outside the space 
of closed loops. (We make no claims of novelty in this 
regard. Various approaches to overco me the sampling 
difficulty have been presented elsewhere! 35 * 36 * 91 ^) 



D. Fluctuating sublattice assignment 

The discussion in the previous section was specific to 
the case in which (i) the AB pattern is regular and (ii) 
the r vectors that have nonzero h(r) only connect sites 
in opposite sublattices. If those conditions hold, there 
are only two possible consequences to the motion of the 
open string: a loop is joined to the string (SNi = — 1) 
or a loop is split off from it (SNi = +1). In both cases, 
represented in Fig.|3]by panels (a)— >(b) and (d)— >(e), the 
AB pattern itself is left undisturbed. 

More generally, as the open string propagates it lays 
down a chain of singlet bonds whose alternating site la- 
bels may be at odds with the traversed sites' current AB 
assignments. A simple workaround is to flip the sub- 
lattice labels as required to correct the mismatch. The 
relevant processes are now those in which a moving open 
string absorbs a closed loop (SNi = —1) or reorganizes it- 
self without impinging on any additional sites (SNi = 0). 
The first case is depicted in Fig. [3] by panels (a)— >(c) and 
the second by (d)— >(f) or (d)— >(g). A crucial consider- 
ation is that, since the singlets are directional, flipping 
sublattice labels along a loop segment has the effect of 
reversing a chain of singlet bonds. If an odd number of 



singlets is effected, the overall sign of the wavefunction 
will change. This is true for all SNi — worm steps. 

The sublattice mismatch can either be a temporary 
condition — lasting only until the worm updates succeed 
in laying down a global AB pattern that is an invari- 
ant of the worm motion — or it may be that the mo- 
tion described by a given h(r) is incompatible with any 
static AB site labelling. For example, consider the one- 
parameter family of short range states on the square 
lattice described by /i(±l,0) = h(0, ±1) = cos 9 and 
h(±l,±l) = smO (with < 9 < tt/4). Regardless of 
the initial sublattice labelling — it can be any random as- 
signment having an equal number of A and B labels — 
the simulation will dynamically establish the checker- 
board pattern provided that 9 — 0. We keep track 
of the AB labelling pattern by measuring a function 
A(Q) = E r ,r< e 4Q -( r ~ r ')(A(r)A(r')), where A(r) takes the 
value —1 or 1 depending on the current sublattice as- 
signment at site r. If 9 = 0, A(Q) starts off broad but 
systematically flows toward the distribution with a single 
delta function peak at Q = (tt, n); once that is achieved, 
the pattern ceases to evolve. Similar behaviour is exhib- 
ited at 9 = tt/4, where the system settles into a static 
pattern with either Q = (n, 0) or Q = (0, 7r). Only in 
those two extreme cases is the sublattice pattern eventu- 
ally static and the simulation sign-problem free. 



III. RESULTS 

As a test of the worm implementation, we compare its 
output to analytical results obtained for the 4x4 lattice. 
We exploit the fact that the bipartite valence bond basis 
Vab for 2N spins is isomorphic to the set of permutations 
on N elements.^ Hence, the basis states have a natural 
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FIG. 4. (color online) (a),(b) Spin correlations Ci and Ci 
between nearest- and next-nearest-neighbour spins, computed 
for the RVB trial wavefunction as a function of the amplitude 
ratio x — h(2, l)//i(l,0). As x becomes increasingly negative, 
numerical evaluation becomes dominated by noise from the 
sign problem, (c) The energy-optimized value x op t remains 
positive up to a g — J2/J1 = 0.40756. (d) The optimized 
trial state gives a good approximation to the true ground 
state energy up to where the Marshall sign rule breaks down. 



lexical ordering via the Lehmer cod d 87 * 88 -* and can easily 
be enumerated. For 4 x 4 = 16 sites, the total number of 
the states is only 8! = 40320, which means that expec- 
tation values of the trial wavefunction can be evaluated 
exactly at very little computational cost. Moreover, we 
can carry out the calculation symbolically. Each observ- 
able takes the form of a rational function of order [16,16]: 



(G) = 



Q(g) 

Z{x) 



o k x" 



T,i=o z i x 



(13) 



The argument of the polynomials appearing in the nu- 
merator and denominator is the real-valued ratio x = 
h(2,l)/h(l,0), and the coefficients are all integers. Co- 
efficients for various observables are listed in Table U 

For this test we have focussed on the nearest- and 
next-nearest-neighbour spin correlation functions C\ — 
J2(i,n\{ s i ■ s j) and c 2 = £«i,,-»<Si • S 3 -); the Q = 



(it, tt) staggered and Q = (it, 0) stripe magnetization 
via M 2 = ^/(-l)^' 1 '"' °<Sr • S r ,>; and D% = 

^Er, r '(-l) ex - (r+r,) ((S r ■ S r+e J(S r , • S^+eJ), the Or- 

der parameter for a columnar dimer crystal. We have 
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FIG. 5. (color online) The level crossings are plotted versus 
1/L 3 and extrapolated to the thermodynamic limit. Several 
different second-order polynomial fits (two shown) are used to 
estimate the uncertainty in the intercept. The solid, blue line 
(fit 1) is an attempt to fit the L > 6 data to Co exp(cia;-r-C23; 2v 
the dashed, green line (fit 2) is a fit to Co + C\L~ i + ciL~ 
L > 4. Our analysis suggests a value g C 2 = 0.5891(3). 
upper inset shows the analysis behind the L — 16 data point, 
which is marked in the main graph as an open circle. The 
lower inset is a magnification of the shaded region. 



verified that the the worm algorithm, conventional bond 
swap Monte Carlo, and exact evaluation give consistent 
results for all these quantities. The comparison for the 
energetics is shown in Fig. |4j Note that in Figs. Qa) 
and Qb) , the stochastic evaluation continues to work in 
some range of x < but quickly breaks down as x be- 
comes strongly negative. The optimized value x op t = 
is achieved for the coupling J2/J1 — 0.40756, and we un- 
derstand this to signal where the Marshall sign rule first 
fails. 

Having established confidence is our numerical imple- 
mentation, we proceed with unbiased optimization cal- 
culations using a static sublattice assignment on lattices 
up to size L — 32. Convergence is limited by statisti- 
cal uncertainty in the (energy to bond count) correlation 
function that determines the local energy gradient,^ and 
it is difficult to reliably optimize for larger system sizes. 
We first consider the checkerboard AB pattern. At g = 0, 
the bond amplitudes are given an initial value of 



h(x, y) = [min(:r, L — a;) 2 + min(y, L — y) 



-3/2 



(14) 



for |x| + |y| odd and zero otherwise. The new set of ampli- 
tudes obtained from this first run serves as the input for 
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FIG. 6. (color online) (Upper panel) The energy per site ver- 
sus the coupling ratio. The solid lines are best energies from 
the trial wavefunction optimization. From bottom to top, sys- 
tem sizes L = 4,6, 8, 10, 12, 16, 20, 24, 28, 32 are shown. The 
solid black line, providing an upper envelope to the curves, 
represents the extrapolation to L — oo. The energies of 
the short range RVB state and of the Heisenberg model on 
the one-dimensional chain are shown for comparison. (Lower 
panel) Magnetization data is shown, with the same system 
sizes now increasing from top to bottom. The solid black lines 
above the grey shading are the L — oo extrapolation. The 
magnetic order for Q = (n, it) and Q = (ir, 0) both vanish in 
the small region between g c i = 0.54(1) and g C 2 = 0.5891(3). 



the next optimization process. That is to say, we daisy 
chain the calculations, at each step using g to seed the 
simulation at g + Sg. An analogous procedure is carried 
out for the stripe AB pattern, starting from g — oo and 
stepping the relative coupling down. 

One finds that the two sets of simulations do not join 
smoothly but instead meet with strongly opposite slopes 
d(H) op t/dg. A detailed extrapolation to the thermody- 
namic limit, presented in Fig. [5j puts the location of the 
energy level crossing at g C 2 = 0.5891(3). As Fig. [6]makes 
clear, this point represents the rightmost edge of an in- 
termediate phase that is magnetically disordered. The 
leftmost edge sits at g c \ = 0.54(1), where the Q = (tt, 7r) 
antiferromagnetism vanishes in a continuous fashion. 

An important detail is that the optimizations are car- 
ried out with the bond amplitudes constrained to have x- 



and y-axis reflection symmetry but not 90° rotation sym- 
metry. In the case of the checkerboard simulation, the 
amplitudes nonetheless realize the full lattice symmetry 
under optimization up to large values of the relative cou- 
pling. For small lattice sizes L — 4, 6, 8, the symmetry 
breaks down beyond values g ~ 0.51,0.55,0.57. For all 
larger sizes, that point is pushed well to the right of g C 2- 
This means that h(r) shares a common symmetry across 
both the staggered magnetic phase and the disordered 
intermediate phase. But it experiences a sudden break 
at the onset of stripe magnetic order, dropping from C4 
to C 2 . 

In the vicinity of g = 0, the optimized bond ampli- 
tudes are positive definite and an almost perfect func- 
tion of bond length. As the frustration increases, the 
amplitudes begin to deviate from circular symmetry, de- 
veloping strong lobes of weight along the x and y axes. 
Bonds not aligned along those preferred directions be- 
come increasingly short ranged, and the eight knight's 
move bonds, equivalent to h(2, 1), start to trend through 
zero to negative values. The extrapolation shown in 
Fig. [7] pinpoints the breakdown of the Marshall sign rule 
at <7mi = 0.398(4). What this suggests is that there is 
strict adherence to a checkerboard Marshall sign rule only 
below (7m 1; in the range <7mi < 9 < 5c2, the sign rule is 
violated even though the overall sign structure is still 
partially consistent with the checkerboard pattern. Wc 
find that the behaviour on the large coupling side is not 
comparable. There, the coupling at which bonds first go 
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FIG. 7. (color online) The coupling strength at which 
the amplitude h(l,2) extrapolates to zero is plotted against 
the inverse linear system size. The shaded region represents 
the envelope containing plausible fits. We estimate that the 
checkboard Marshall sign rule fails at jmi = 0.398(4) in the 
thermodynamic limit. 
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in each case starting from a random AB pattern and a 
random loop configuration. The bond amplitudes are 
initialized with h(r) forming a broad peak around r = 
and having no zero entries. We perform a crude simula- 
tion in which any sign associated with the worm moves 
is thrown away. Otherwise, the optimization of h(r) pro- 
ceeds as before. What we find is a result that exactly 
tracks the state of lower energy produced by assuming 
one of the two static AB patterns. The simulation flows 
to the checkboard for all g < g C 2 and to the stripe for all 
9 > 9c2] the peak in A(Q) jumps discontinuously. Ob- 
viously we should not read too much into a result that 
follows from an uncontrolled approximation (sampling by 
ignoring the signs), but it does give us a sense that the 
stability of the checkerboard pattern through the inter- 
mediate phase and the abrupt change in Marshall sign 
structure at g C 2 might be genuine features of the model. 

Given our choice of trial wavefunction, the intermedi- 
ate phase that we compute is certainly not a bond crys- 
tal. The dimer correlations are somewhat enhanced in 
the strongly frustrated region, but with increasing lattice 
size they show clear convergence to zero. Still, spatially 
resolved dimer correlations do give us important infor- 
mation. One can see in Fig. [9] that the optimized state 
shows the same pattern of dimer correlation and anticor- 
relation as the short-bond-only RVB but it decays much 
faster as a function of dimer seperation. The comparison 
is made more explicit in Fig. |10[ which shows correla- 
tions along a line and a stack of dimers. The functions 
measured are 



FIG. 8. (color online) Schematic representation of the 
model's zero temperature phase diagram. Critical couplings 
g c i and g C 2 mark the boundaries of the magnetically disor- 
dered phase. Staggered order ends with a continuous transi- 
tion at <? c i; stripe order ends with a first-order transition at 
g C 2- The three diagrams on the right illustrate the optimized 
h(r) values at g = 0, g — 0.55, and g = 0.8. Each circle, offset 
by a vector r measured from the centre, has an area propor- 
tional to the corresponding h(r) value. The tags on the left 
indicate the Marshall sign structure that predominates. 



negative scales as <?M2 ~ 1/L 4 and hence does not con- 
verge in the thermodynamic limit. We interpret this to 
mean that the static stripe pattern is only ever a weak 
description of the Marshall sign structure. See Fig. [8] 

We have attempted to confirm this picture by running 
simulations in which the Marshall sign structure is de- 
termined dynamically. More specifically, we want to ver- 
ify that the strongly first-order transition at g C 2 is not 
merely an artifact of two incompatible sublattice conven- 
tions colliding. If permitted, might the system smoothly 
interpolate over some range of <?, with the peak in A(Q) 
migrating from (ir, it) to (ir, 0)? We follow the procedure 
outlined in Sect. II D[ whereby the sublattice labelling is 
no longer fixed and the worm motion itself is allowed to 
reconfigure the current AB pattern. Our approach is to 
simulate for various g values — with no daisy chaining — 



Ciinc(rf) = (D(0,0)D(d,0)) - (D(0,0))(D(d,0)) 
C sta ck(d) - (D(0,0)D(0,d)) - (D(O,O))0(O,d)>, 



(15) 



which we have expressed in terms of the x-directed bond 
operator D{x, y) = S(x, y) ■ S(x + 1, y). 



IV. 



CONCLUSIONS 



We have used an optimized valence bond trial wave- 
function to study the square-lattice J1-J2 Heisenberg 
model, with an eye to both mapping out the zero- 
temperature phase diagram and determining how the 
Marshall sign structure breaks down near the phase 
boundaries. In the first instance, we fix the AB sublattice 
labelling to coincide with the order that exists at small 
and large coupling. For each lattice size, the intermedi- 
ate phase is approached in two independent simulations 
(or, rather, chains of history-dependent simulations) by 
evolving the states progressively out of the two ordered 
phases, minimizing their energy at each step. 

Finite-size scaling of the dimer order parameter sug- 
gests that there is no long-range dimer order for any 
range of g values. This is to be expected, since the trial 
state explicitly ignores bond-bond correlations. Measure- 
ments of the staggered magnetization show clear evidence 
of a continuous phase transition in which the staggered 
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tional symmetry. (Since the trial state is least able to de- 
scribe the intermediate phase — again, because of its lack 
of bond-bond correlations — we should probably view g c i 
and <?c2 as upper and lower bounds, respectively, on the 
true positions of the phase boundaries.) We have also 
performed calculations in which no sublattice labelling is 
put in by hand. Our results suggest that the checker- 
board AB pattern is the correct choice throughout the 
intermediate phase. Hence, within the context our trial 
wavefunction scheme, we surmise that the state beyond 
g c i is a "bosonic" gapped spin liquid, and the gap to its 
excitations has a minimum at (n,n). 

Figure [8] gives a quick summary of our results. We ob- 
serve that at high frustration the bond amplitudes take 
on a highly anisotropic form. This is quite different from 
the long-bond to short-bond picture that is usually in- 
voked. Recall that Liang, Ducot, and Anderson studied 
long range RVB states on the square lattice with ampli- 
tudes h ~ r~ p that decay as a powerlaw in the bond 
length rP^ In that framework, the state becomes mag- 
net ically d isordered when p exceeds a critical value of 
3 3j 30|85|8 6] anc j ^- ne en y re family of states in the range 

p > 3.3 is continuously connected to p = 00, which is the 
short-bond-only RVB. The intermediate phase state ob- 
tained in our simulations is of a quite different character. 
The state is magnetically disordered not because its bond 
amplitudes are sufficiently short ranged but because they 
are sufficiently anisotropic. 

This work was supported by a Discovery grant from 
NSERC of Canada. 



FIG. 9. Grid lines depict the dimer correlations Cijki = 
{(Si ■ Sj)(Sfc • Si)) on the nearest-neighbour links (k, I) of the 
square lattice, measured with respect to the thick, dark dimer 
(i,j) at the centre. The correlations are computed for the 
L — 28 system. The grayscale intensity represents correlation 
strength [presented as the fourth power of (1 + 1.5r )(7], 
and dotted lines indicate a negative (anticorrelated) value. 
The top panel shows results for the purely short ranged RVB 
state, presented for comparison's sake. The bottom panel 
shows results for the energy-optimized state at g = 0.58. In 
each 10 x 10 piece of the of the full valence bond loop 

configuration, obtained from a snapshot of the Monte Carlo 
simulation, is overlaid. 
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magnetization vanishes at g c \ = 0.54(1). On the other 
edge of the intermediate phase, an energy level cross- 
ing at g C 2 = 0.5891(3) results in the sudden disappear- 
ance of the otherwise robust stripe magnetization. This 
is accompanied by the restoration of the system's rota- 



FIG. 10. (color online) Dimer correlations of the product- 
amplitude trial wavefunction optimized at g = 0.58 (solid 
points) and the short-bond-only RVB state (open points) are 
compared on the L — 28 lattice. Presented are the dimer line 
(squares) and dimer stack (circles) correlation functions. See 
the text for definitions. 
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n 


Ci 


c 2 


Z 


M 2 (-K,n) 


M 2 (tt,0) 


D 





22241280 


9902080 


1559232 


113983488 


17383424 


4376064 


1 


194568192 


104726528 


13008384 


1117618176 


139902976 


28540928 


2 


997232640 


585695232 


66018816 


6104383488 


709410816 


127591424 


3 


3395051520 


2137292800 


223842816 


21861335040 


2381496320 


389861376 


4 


8564477952 


5689352192 


568694016 


57653526528 


6069354496 


932687872 


5 


16547069952 


11594661888 


1108661760 


116342292480 


11792498688 


1697314816 


6 


25797685248 


18932629504 


1767412224 


189239033856 


18888998912 


2580870144 


7 


32679444480 


25148850176 


2302253568 


250229981184 


24519589888 


3165620224 


8 


34418749440 


27661209600 


2528419968 


275349995520 


27030159360 


3329164288 


9 


29878050816 


25148850176 


2302253568 


250229981184 


24519589888 


2857185280 


10 


21512073216 


18932629504 


1767412224 


189239033856 


18888998912 


2089987072 


11 


12538503168 


11594661888 


1108661760 


116342292480 


11792498688 


1223688192 


12 


5848903680 


5689352192 


568694016 


57653526528 


6069354496 


594391040 


13 


2070282240 


2137292800 


223842816 


21861335040 


2381496320 


218601472 


14 


528863232 


585695232 


66018816 


6104383488 


709410816 


60980224 


15 


84836352 


104726528 


13008384 


1117618176 


139902976 


11331584 


16 


6254592 


9902080 


1559232 


113983488 


17383424 


1074688 



TABLE I. The integer coefficients for the rational polynomial that expresses various measured quantities for the product 
amplitude trial state on the 4x4 lattice. The columns correspond to the nearest- and next-nearest-neighbour spin correlations, 
the partition function (denominator), the staggered and stripe magnetization, and the columnar dimer order parameter. 



